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Complex  Variable  Methods  for  Shape  Sensitivity  of  Finite  Element  Models 

Andrew  Voorhees,  Ronald  Bagley,  Harry  Millwater 


Abstract 

Complex  variable  methods  have  some  potential  advantages  over  classical  finite 
differencing  methods  for  sensitivity  analysis.  Two  methods,  complex  Taylor  series  expansion 
and  Fourier  differentiation,  are  applied  and  compared  to  central  differencing  for  shape  sensitivity 
analysis.  A  two  dimensional  finite  element  model  with  an  analytical  solution  is  chosen  for  use  in 
comparing  the  accuracy  of  the  methods.  It  is  found  that  for  the  accuracy  of  the  model  chosen, 
the  error  in  the  sensitivities  is  primarily  defined  by  the  error  in  the  solution,  not  the  error  in  the 
sensitivity  method. 

Introduction 

Finite  element  analysis  is  a  powerful  tool  in  computational  engineering.  The  method 
allows  for  the  approximate  solution  of  boundary  value  partial  differential  equations.  Given 
enough  computational  power  it  is  possible  to  solve  highly  non-linear  equations,  on  extremely 
complicated  and  intricate  domains,  with  arbitrary  boundary  conditions  [1].  Finite  element 
analysis  has  become  an  important  part  of  the  design  process  for  mechanical  engineers.  Accurate 
numerical  solutions  can  reduce  the  need  for  costly  laboratory  experiments.  Furthermore,  the 
effect  that  a  small  design  change  has  on  the  performance  of  the  model  can  be  calculated  quickly, 
avoiding  the  need  to  build  new  test  specimens. 

The  effect  that  a  small  design  change  has  on  the  output  parameters  of  the  model  can  be 
characterized  through  a  process  known  as  sensitivity  analysis.  A  sensitivity  is  in  fact  identical  to 
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a  partial  derivative,  in  that  it  quantifies  the  change  in  one  parameter  due  to  a  change  in  another 
parameter,  with  all  other  inputs  held  constant.  This  means  that  sensitivity  analyses  can  be 
performed  using  simple  numerical  differentiation  techniques.  For  design  work,  the  most 
important  type  of  sensitivity  is  the  shape  sensitivity  [2].  A  shape  sensitivity  tells  the  designer 
what  effect  a  small  change  in  the  size  or  shape  of  the  domain  will  have  on  the  output.  Examples 
of  shape  sensitivities  include  the  sensitivity  of  the  lift  generated  by  an  aircraft  wing  to  the  shape 
of  the  wing’s  cross  section,  or  the  sensitivity  of  the  stress  in  a  beam  due  to  a  change  in  it’s 
length. 

Sensitivities  are  easily  calculated  through  numerical  differentiation  techniques. 
Traditionally,  the  numerical  differentiation  method  of  choice  has  been  finite  differencing.  Finite 
differencing  requires  that  a  function  be  evaluated  at  additional  sample  points  and  the  derivative 
of  the  function  can  be  evaluated  by  calculating  the  difference  in  the  function’s  value  between 
two  of  the  sample  points.  Central  differencing  (CD)  is  a  widely  used  form  of  finite  differencing 
chosen  for  it’s  increased  accuracy.  Over  the  last  ten  years,  alternative  numerical  differentiation 
techniques  have  emerged  for  use  in  sensitivity  analysis.  Two  of  these  methods  are  complex 
Taylor  series  expansion  (CTSE),  also  referred  to  as  the  complex  step  derivative  method,  and 
Fourier  differentiation.  These  methods  offer  more  accurate  and  stable  derivatives  compared  to 
CD. 

CTSE  was  first  described  by  Lyness  and  Moler  in  the  late  1960’s  [3,4].  It  reemerged  as  a 
tool  for  engineering  analysis  with  a  paper  by  Squire  and  Trapp  in  1998  [5].  Since  then  it  has 
been  used  in  a  wide  variety  of  engineering  fields  including  computational  fluid  dynamics, 
dynamic  system  optimization  and  many  more  [6-13].  In  all  of  these  fields  CTSE  has  offered  a 
great  improvement  in  accuracy  over  CD.  However,  CTSE  was  only  found  to  offer  similar 
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accuracy  to  CD  for  use  in  the  calculation  of  shape  sensitivity  problems  for  one  dimensional  and 
two  dimensional  finite  element  models  [13]. 

FD  was  also  developed  by  Lyness  in  the  late  60’s  and  early  70’s[3,4, 1 5] .  The  method 
has  was  further  described  by  Henrici  and  more  recently  Bagley  [16,17].  The  method  utilizes 
additional  sample  points  in  the  complex  plane  and  an  FFT  routine  to  calculate  derivatives 
including  high  order  derivatives  with  exceptional  accuracy.  To  date  the  method  has  not  been 
widely  used  for  the  determination  of  sensitivities  for  engineering  problems. 

Although  the  scope  of  this  paper  does  not  cover  them,  several  other  sensitivity  methods 
that  do  not  rely  on  the  evaluation  of  sample  points  have  come  into  use  in  the  engineering 
community  [18-24].  These  codes  are  typically  more  difficult  to  apply  to  problems  in  that  they 
require  extra  coding  or  the  derivation  and  solution  of  new  equations  rather  than  simply 
generating  extra  sample  points  and  applying  simple  numerical  differentiation  formulas.  For  the 
most  part  they  are  restricted  to  problems  in  which  high-dimensional  gradients  need  to  be 
calculated.  The  first  of  these  methods  is  automatic  differentiation.  This  method  is  based  on  the 
concept  of  the  chain  rule,  and  the  fact  that  a  large  code  thought  of  as  a  single  function  composed 
of  several  small  functions  each  having  its  own  partial  derivative  [18].  Thus  by  tracking  the  use 
of  each  small  function  in  a  code,  i.e.  every  multiplication  or  subtraction  operation,  and  storing 
the  derivative  information,  sensitivities  can  be  obtained  by  carefully  applying  the  chain  rule. 
Fortunately,  several  automatic  differentiation  codes  for  a  variety  of  programming  languages 
already  exist.  These  include  ADIFOR  and  ADI-C  codes  for  FORTRAN  and  C  respectively 
[25,26].  This  method  has  been  found  to  be  more  accurate  than  the  sampling  methods  because 
each  derivative  evaluation  can  be  done  symbolically,  which  means  that  the  error  will  due  to 
machine  round-off.  Automatic  differentiation  is  also  very  efficient  for  computing  high- 
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dimensional  gradients.  Direct  differentiation  and  the  adjoint  method  are  two  additional  methods 
[20].  In  direct  differentiation,  sensitivities  are  calculated  through  linearization  and  discretization 
of  the  original  differential  equations  to  produce  discrete  linearized  state  equations  that  can  be 
solved  to  produce  the  desired  sensitivities.  The  adjoint  method  is  similar,  except  the 
linearization  is  replaced  with  an  adjoint  operation.  These  methods  have  been  found  to  be  very 
useful  in  the  field  of  CFD  where  problems  often  involve  large  numbers  of  input  variables  [20] . 

The  goal  of  this  paper  is  to  demonstrate  the  use  of  complex  variable  methods  for  the 
calculation  of  shape  sensitivities  for  a  simple  two-dimensional  finite  element  model.  The  two 
dimensional  plane  stress  elastic  model  of  a  thick  walled  cylinder  under  uniform  boundary 
pressure  will  be  used  as  a  numerical  example.  This  problem  has  an  analytical  solution,  which 
can  be  differentiated  to  determine  the  shape  sensitivities.  The  sensitivities  calculated  by  CTSE 
and  FD  will  be  compared  to  the  analytical  sensitivities  as  well  as  those  calculated  by  CD. 


Methodology 


Numerical  Differentiation 

Numerical  differentiation  is  a  process  through  which  an  estimate  of  a  function’s 
derivative  can  be  obtained.  A  derivative  is  defined  as  the  limit  of  the  change  in  a  function’s 
value  across  two  different  points,  as  the  distance  between  the  two  points  goes  to  zero. 


f\xo)  =  lim 

X~>X0 


f(x)~  fix  J 
x-xa 


(1) 


Finite  differencing  methods  calculate  derivatives  by  estimating  the  limit  in  eq.  1  as  a  difference 
between  a  function  evaluated  a  two  distinct  points  located  a  distance  h  apart. 

fixo  +  h)-fixo ) 


f'ixj: 


h 


(2) 


4 


This  distance,  h,  is  often  called  the  step  size.  When  h  is  positive,  the  method  is  referred  to  as 


forward  differencing.  When  h  is  negative  it  is  called  backwards  differencing.  When  the  forward 
difference  and  the  backwards  difference  are  averaged,  the  method  is  called  central  differencing. 
The  equation  for  central  differencing  is  as  follows. 

fuY  \  ^.f(x0  +  h)-f(x0-h) 

o)~  2/7  (3) 

The  approximation  of  the  derivatives  as  a  difference  between  two  non-identical  numbers  leads  to 
error  due  to  the  truncation  of  terms  in  the  functions  Taylor  series.  This  error  can  be  eliminated 
by  making  the  step  size  as  small  as  possible.  However,  as  the  step  size  gets  very  small,  a  new 
source  of  error  arises.  This  new  error  is  round-off  error  and  it  is  due  to  the  fact  that  a  computer 
cannot  accurately  calculate  a  small  difference  between  two  large  numbers.  This  means  that  for 
finite  differencing  there  is  a  lower  limit  on  the  step  size  and  also  a  limit  on  the  maximum 
achievable  accuracy. 

For  the  forward  differencing  method  all  Taylor  series  terms  above  the  first  order  term  are 
ignored.  This  means  that  the  order  of  accuracy  for  a  given  step  size  is  0(h).  By  using  CD  all  the 
even  order  terms  in  the  Taylor  series  cancel  out  and  the  accuracy  of  the  method  becomes  0(h  ). 
The  increased  accuracy  of  CD  is  the  reason  it  has  become  the  standard  method  for  sensitivity 
calculations. 

Higher  order  derivatives  can  also  be  calculated  through  CD,  by  using  additional  sample 
points.  The  formula  for  the  second  derivative  is. 

f(x0  +  h)-2f(xa)  +  f(xa  -  h) 


f{2)M 


Jr 


(4) 


Formulae  exist  for  derivatives  above  second  order  but  are  not  printed  here.  One  of  the  problems 
with  CD  is  that  the  calculation  of  higher  order  derivatives  requires  more  sample  points  and  more 
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difference  operations.  Each  additional  difference  operation  results  in  an  increase  in  the  round-off 
error,  which  further  restricts  the  lower  limit  of  h.  This  means  that  CD  is  not  a  good  choice  for 
the  calculation  of  higher  order  derivatives. 

CTSE  is  another  numerical  differentiation  method.  CTSE  uses  the  orthogonality  of  the 
real  and  imaginary  axes  of  the  complex  plane  to  calculate  derivatives  with  fewer  difference 
operations  and  in  turn  less  round-off  error  when  compared  to  CD.  CTSE  is  capable  of 
calculating  both  the  first  and  second  order  derivatives  from  a  single  sample  point  located  at 
x0+ih.  The  formulae  for  the  derivatives  can  be  derived  from  the  Taylor  series  representation  of 
the  function  evaluated  at  the  complex  sample  point. 


f(x„  +  ih)  =  f(xa)  +  fl\xa)  ■  l-^~  +  /(2,(x0)  • 


+  /(3)(*> 


('■hf 

3! 


+ . 


(5) 


Taking  the  imaginary  part  of  both  sides  of  eq.  5  and  solving  for  the  first  derivative  will  result  in 
an  approximation  with  accuracy  Oih2). 

/W)Cm(/'fc  +  ft»  (6) 

n 


It  is  noted  that  for  the  first  derivative  no  difference  operation  is  needed.  This  means  that  the  step 
size  can  be  made  arbitrarily  small  with  no  concern  about  increasing  round-off  error.  Taking  the 
real  part  of  eq.  5,  the  formula  for  the  second  derivative  with  error  0{h  )  can  be  derived. 

f(x)  x  2(f(Xo)-Mf(xo  +  ih)))  (7) 


It  is  noted  that  the  second  derivative  contains  a  difference  operation  meaning  that  round-off  error 
will  be  a  problem  if  h  is  set  too  small.  By  using  more  sample  points  it  is  possible  to  solve 
equation  5  to  obtain  the  higher  order  derivatives. 

FD  can  be  derived  from  the  Cauchy  integral  formula. 
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(8) 


f(n\z) 


2m v 


m 

(^-z)”+1 
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which,  relates  any  nth  order  derivative  to  a  simple  contour  integral  in  the  complex  plane.  The 
form  of  this  integral  is  identical  to  the  form  of  a  Fourier  integral,  and  hence  it  can  be  evaluated 
using  an  FFT  routine.  If  the  function  of  interest  is  evaluated  at  N  sample  points  along  a  circular 
contour  in  the  complex  plane  centered  on  the  initial  point,  a  vector  of  the  sampled  data  can  be 
run  through  an  FFT  routine,  and  the  output  will  be  the  first  N  terms  in  the  functions  Taylor 
series.  The  nth  order  derivative  of  the  function  can  then  be  calculated  from  the  Taylor  series 
coefficients  by  using  the  following  relationship. 

ay. 

h"  (9) 

For  more  information  on  Fourier  differentiation  see 

Bagley,  2006  [17]. 


J  \  os 


where  a„  is  the  nth  Taylor  series  coefficient. 


Finite  Element  Analysis 

A  two-dimensional  finite  element  code  capable  of  solving  the  2-D  equations  of  elasticity 
under  the  assumptions  of  plane  stress  was  written  in  Matlab.  The  code  uses  second  order  shape 
functions  and  six -node  triangular  elements.  The  equations  for  the  shape  functions  in  terms  of  the 
natural  coordinates  can  be  found  in  Huebner,  on  page  543  [27].  The  use  of  natural  coordinates 
allows  for  all  numerical  integration  to  be  replaced  with  algebraic  manipulation,  which  greatly 
improves  the  run  time  of  the  program.  Geometry  and  meshes  were  created  using  the  Comsol 
finite  element  package  and  were  then  imported  into  Matlab.  After  importation,  the  edges  of  the 
elements  were  reconstructed  so  that  all  element  edges  were  linear,  and  all  non-vertex  nodes  were 
placed  on  the  midpoint  of  the  edges.  The  code  used  a  conjugate  gradient  solver  with  a  minimum 
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error  of  le-8.  After  solving  the  global  system  of  equations  for  the  displacements,  the  stresses  at 
each  of  the  nodes  were  calculated  on  an  element  by  element  basis  using  the  equations  of 
elasticity,  second  order  shape  functions  and  the  nodal  displacements.  Since,  the  nodal  stress 
calculated  from  one  element  does  not  always  agree  with  the  nodal  stress  calculated  using  another 
element,  stress  averaging  was  used  to  smooth  the  stress  field.  The  code  then  calculated  the 
principal  stresses  and  the  Von  Mises  stresses  using  the  stresses  at  the  nodal  points. 

In  order  to  perform  the  shape  sensitivity  calculations,  a  small  step  must  be  added  to  the 
nodal  coordinates.  This  was  done  by  first  identifying  every  node  that  was  located  on  the 
geometry  feature  that  was  being  sampled.  For  instance,  if  the  sensitivity  of  the  solution  to  the 
radius  of  a  hole  was  being  calculated,  the  program  would  identify  each  node  that  was  located  on 
the  surface  of  the  hole.  The  code  would  then  add  a  small  step  (complex  or  real  depending  on  the 
method)  to  the  nodal  coordinate  of  each  of  the  identified  nodes.  All  other  nodes  in  these 
elements  would  remain  the  same.  If  the  step  size  chosen  for  the  numerical  differentiation 
technique  is  very  large,  be  it  real  or  imaginary,  the  elements  may  become  distorted,  which  can 
lead  to  poor  numerical  results,  and  as  such  the  step  sizes  and  sampling  radii  of  the  differentiation 
methods  were  typically  kept  much  smaller  than  the  edge  lengths  of  the  elements. 

Numerical  Example:  Thick  Walled  Cylinder 

In  order  to  test  the  accuracy  of  the  three  numerical  differentiation  techniques  it  is 
necessary  to  find  a  problem  with  an  analytical  solution.  The  thick  walled  cylinder  under  uniform 
boundary  pressure  is  such  a  problem.  The  equations  that  govern  the  stress  through  the  thickness 
of  the  cylinder  are  given  in  eq.  10  [28], 


(10) 


2  2  i  2 

_  _  n  f2  P  1  r2  p 

ur  222  22 

r2  -n  r  >  2  -n 

2  2,  2 

_ n  r2  p  1  r2  p 


1  2 

°e=-^T- 


2  2  2  2  2 
r2  ,  r  r2  -rx 


where,  r\  is  the  inner  radius,  o  is  the  outer  radius,  and  p  is  the  boundary  pressure.  For  this 
example,  the  inner  radius  of  the  cylinder  is  0.5  m,  the  outer  radius  is  1  m  and  the  boundary 
pressure  is  10  kPa.  The  stress  equations  given  in  eq.  10  can  be  differentiated  with  respect  to  the 
inner  radius  to  generate  the  sensitivities  of  the  stresses.  The  sensitivities  of  the  stresses  with 
respect  to  the  inner  radius  appear  in  eq.  1 1  and  the  second  order  sensitivities  appear  in  eq.  12. 


for  _  ry  >V22P  \  r2  1 

o  Z  ,  2  2,.2  2  1 
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Vi  P  r2 
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The  problem  was  solved  using  four  different  meshes  in  order  to  examine  the  convergence 
of  the  error  in  the  solution.  The  coarsest  mesh  contained  888  elements,  and  1868  nodes,  the  next 
mesh  contained  1792  elements  and  3716  nodes,  the  third  mesh  contained  6184  elements  and 
12,608  nodes,  and  the  finest  mesh  consisted  of  25,124  elements  and  50,724  nodes.  The  solutions 
for  the  radial  and  tangential  stresses  as  calculated  using  the  6184  element  mesh  are  shown  in 
figure  1.  The  following  norm  was  selected  in  order  to  compare  the  error  in  the  four  different 


mesh  cases. 
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error  =  ■ 


mean(  cr 


analytical  ^ numerica 


l) 


rneani  (Janalytlca\) 


(13) 


The  following  formula  was  used  to  calculate  the  error  plotted  in  the  figures  for  this  example. 


error = 


analytical  numerical 


maX(  ^analytical) 


(14) 


This  formula  was  chosen  so  as  to  minimize  the  influence  of  large  errors  at  locations  where  the 
solution  is  near  zero,  such  as  at  the  inner  radius  for  the  radial  stress.  Table  1  shows  the  norm  of 
the  error  in  both  the  radial  and  tangential  stress  solutions  for  each  mesh  case.  This  data  is  shown 
graphically  in  figure  2.  It  is  seen  that  each  successive  mesh  iteration  reduces  the  error  by 
approximately  half  an  order  of  magnitude.  The  amount  of  computational  time  (wall  time) 
needed  to  solve  each  mesh  case  is  shown  in  table  2.  Given  the  amount  of  computational  time 
required  for  25,124  element  case  and  the  fact  that  the  complex  sensitivity  solutions  will  require 
three  times  more  computation  than  the  real  valued  case,  the  6184  element  case  was  used  to 
generate  the  first,  second  and  third  order  sensitivities.  For  each  sensitivity  method  the  step  size 
or  sampling  radius  was  0.001  which  is  approximately  l/30th  of  the  average  element  edge  length. 
CTSE  and  CD  were  both  performed  using  as  few  sample  points  as  possible,  and  FD  was 
performed  using  6  sample  points.  The  norm  of  the  error  in  the  sensitivities  appears  in  table  3. 
The  error  in  the  first  order  sensitivities  of  the  radial  stress  over  the  entire  domain  appear  in  figure 
3,  and  the  error  in  the  second  order  sensitivities  of  the  radial  stress  appear  in  figure  4.  These 
figures  show  only  very  slight  differences  between  the  three  methods.  It  is  also  seen  that  along 
the  inner  circumference  of  the  cylinder  the  error  is  very  large.  This  is  due  to  the  fact  that  the 
sensitivity  cannot  be  accurately  calculated  on  the  sampled  surface  itself,  because  the  boundary 
conditions  require  the  solution  to  be  fixed  at  the  inner  surface. 
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The  norm  of  the  error  in  each  case  is  mostly  independent  of  the  method  selected.  This  is 
especially  true  for  the  first  order  sensitivities.  The  third  order  sensitivities  show  some  small 
differences  between  the  methods,  with  CTSE  having  the  highest  accuracy.  The  fact  that  the  error 
is  similar  between  all  three  methods  points  to  the  fact  that  the  error  in  the  solution  is  dominating 
any  errors  arising  from  the  differentiation  methods  themselves.  This  is  seen  by  looking  at  the 
first  and  second  order  sensitivities  of  the  radial  stress  as  a  function  of  the  number  of  elements 
(table  4).  It  is  quickly  seen  that  each  additional  mesh  refinement  increases  the  accuracy  of  the 
method.  This  reduction  in  the  errors  of  the  sensitivities  is  similar  to  the  reduction  in  the  error  of 
the  solution  due  to  further  mesh  refinement  seen  in  table  1. 

The  error  in  the  first  and  second  order  sensitivities  of  the  radial  stress  calculated  using 
three  different  step  sizes  are  shown  in  table  5.  It  is  seen  that  changing  the  step  size  does  not  have 
much  effect  on  the  accuracy  of  the  sensitivity.  This  is  a  further  indicator  that  the  accuracy  of  the 
solution  is  limiting  the  accuracy  of  the  sensitivities,  not  the  accuracy  of  the  numerical 
differentiation  methods.  One  exception  is  the  second  order  sensitivity  at  the  smallest  step  size, 
0.0001  or  l/300th  of  the  average  element  edge  length.  At  this  step  size  each  method  produces 
sensitivities  that  are  less  accurate  than  those  calculated  with  a  larger  step  size.  This  indicates  that 
the  machine  round-off  error  associated  with  this  step  size  may  be  similar  in  magnitude  to  the 
error  due  to  the  solution. 

Numerical  Example:  Disc  in  Diametrical  Compression 

One  of  the  classic  tests  in  material  analysis  is  the  disc  in  diametrical  compression  [20]. 

In  this  test  a  circular  disc  is  loaded  in  compression  along  its  y-axis.  The  load  is  modeled  as  a 
point  load.  This  loading  and  geometry  generates  a  very  nice  uniform  tensile  stress  along  the  x- 
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axis  of  the  specimen.  It  is  thus  useful  in  examining  the  tensile  properties  of  a  material  without 
actually  loading  the  specimen  in  tension.  This  test  is  also  known  as  the  indirect  tensile  test  or  the 
Brazil  nut  test. 


The  diametrical  compression  test  has  an  analytical  solution  that  can  be  derived  through 
simple  superposition.  The  solution  of  the  stresses  is  seen  in  eq.  15  [28]. 


-2  P 

(R-y)x2  |  ( R  +  y)x 2 

n 

_(x2  +  (R-y)2)-  (x2  +  (R  +  y)2j 

-2  P 

C R-y y  ,  (R+yY 

n 

_(x2  +  (R-y)2j-  (x2  +  (R  +  y)2j- 

2  P 

C R-yfx  |  ( R  +  yfx 

71 

(x~  +  (R-T)2)  (x~  +  {R  +  v)2^ 

J_ 

2R 

J_ 

2R 


(15) 


In  these  equations,  P  is  the  magnitude  of  the  point  load,  R  is  the  radius  of  the  disc,  and  x  and  y 
specify  the  location  at  which  the  stress  is  calculated,  with  the  point  (0,0)  located  at  the  center  of 
the  disc.  The  analytic  solutions  of  eq.  15  can  be  differentiated  to  yield  the  sensitivities  with 
respect  to  the  radius  of  the  disc.  The  equations  for  the  first  two  sensitivities  of  the  normal  stress 
in  the  x-direction  with  respect  to  the  radius  are 


d<Jx  _  2  P 


n 


dR 

d2a „  -2 P 


(3 R2  -  6 Ry  -x2  +  3 y2)x2  (3 R2  +  6 Ry  -x2  +  3 y2)x2  1 

~\~  -j  ~Y~ 


( x2  +  (R-y )2) 


(t2  +  (R  +  j)2) 


2  R2 


(16) 


dR 


n 


12(R  - y)(R2  -  2Ry  -  x2  +  y  )x  \2(R  +  y)(R2  +2Ry  —x  2  +  y2)x2 


2\4 


(R  -2 Ry  +  x~  +  y) 


2\4 


(R  +  2 Ry  +  x  +  y  ) 


R 


The  closed-form  solutions  for  the  sensitivities  make  this  problem  another  excellent  choice  for 
exploring  the  use  of  the  complex  variable  sensitivity  methods. 

The  diametrical  compression  test  model  was  solved  using  three  different  meshes,  with  a 
coarse  mesh  consisting  of  1148  elements  and  2357  nodes,  a  moderately  refined  mesh  of  2502 
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elements  and  5093  nodes  and  a  fine  mesh  with  8,374  elements  and  16,909  nodes.  The  solution 
for  the  stresses  as  calculated  using  the  fine  mesh  appears  in  figure  5.  It  should  be  noted  that  for 
each  figure  in  this  example,  no  solution  is  plotted  for  the  elements  that  share  the  node  where  the 
load  is  applied.  This  is  due  to  the  fact  that  the  stress  on  this  node  would  be  infinite.  The  error 
norm  used  for  this  example  is  the  same  as  for  the  thick  walled  cylinder  example.  The  formula  for 
the  error  plotted  in  the  figures  is. 


error = 


^ analytical  ^ numerical 


(J 


analytical 


(17) 


Table  6  shows  the  error  for  the  three  stresses  for  the  three  different  mesh  sizes.  As  before,  it  is 
seen  that  each  successive  mesh  refinement  results  in  a  fairly  large  reduction  in  the  error  norm. 
This  is  seen  visually  in  figure  6.  The  red  areas  are  regions  of  high  error.  As  the  number  of 
elements  increases,  it  is  seen  that  the  total  size  of  the  red  regions  decreases  significantly  as  the 
mesh  is  further  refined.  The  computational  time  required  to  generate  one  solution  appears  in 
table  7.  The  errors  in  the  sensitivities  of  the  normal  stress  in  the  x-direction,  calculated  by  each 
of  the  three  methods  are  plotted  in  figure  7  (first  order)  and  figure  8  (second  order)  for  the  fine 
mesh  case.  It  is  seen  that  there  is  again  high  error  along  the  circumference  of  the  disc  due  to 
boundary  conditions.  It  is  also  noted  that  a  few  lines  of  high  error  form  inside  the  discs.  These 
lines  represent  regions  where  the  analytical  sensitivity  is  zero  or  near  zero.  Since  the  analytical 
sensitivity  appears  in  the  denominator  of  the  error  formula  given  in  eq.  17,  the  error  becomes 
very  large  when  the  analytical  sensitivity  tends  towards  zero.  The  norm  of  the  error  in  the 
sensitivities  are  shown  in  table  8.  It  is  seen  that  the  norm  of  the  error  is  much  larger  in  this  case. 
This  is  in  keeping  with  the  fact  that  the  error  norm  of  the  solutions  themselves  are  much  higher 
than  for  the  first  example.  This  time  it  is  seen  that  there  is  not  a  large  dependence  of  the  error 
norm  on  the  number  of  elements.  It  is,  however,  seen  that  there  is  not  much  difference  between 
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the  methods  themselves.  This  points  to  the  fact  that  the  error  in  the  derivatives  is  not  due  to  the 
truncation  error  of  the  derivative  methods,  since  the  truncation  error  of  CTSE  and  CD  should  be 
of  the  order  0(h2 )  while  the  truncation  error  of  FD  should  be  0(h6).  The  lack  of  dependence  on 
the  choice  of  method  is  further  shown  in  table  9  where  the  norm  of  the  errors  is  shown  for  the 
2502  element  mesh  for  three  different  step  sizes,  or  sampling  radii.  Very  little  variation  in  the 
error  is  seen  as  a  function  of  the  step  size,  which  is  not  the  behavior  of  truncation  error. 

Conclusions 

CTSE,  FD,  and  CD  can  all  be  used  to  calculate  shape  sensitivities.  This  marks  the  first 
time  that  FD  has  been  used  for  calculating  finite  element  shape  sensitivities.  Unfortunately,  For 
2-D  finite  element  problems  using  second  order  shape  functions,  the  error  in  the  model  is  greater 
than  the  error  due  to  the  truncation  errors  associated  with  the  numerical  differentiation  methods. 
This  means  that  the  complex  variable  sensitivity  methods  do  not  offer  extra  accuracy  compared 
to  CD.  If  the  model  were  made  to  be  more  accurate,  such  as  with  higher  order  basis  functions  or 
more  elements,  then  FD  and  CTSE  could  offer  improved  accuracy.  It  has  been  shown  that  FD  is 
capable  of  producing  highly  accurate  sensitivities  for  functions  with  solutions  that  are  accurate  to 
machine  precision  [17].  The  trade-off  for  the  increased  accuracy  of  FD  is  the  requirement  of 
several  complex  sample  points.  It  takes  three  times  more  computational  effort  to  generate  a 
complex  sample  than  a  real  valued  sample.  Thus  for  the  problems  described  in  this  paper,  FD  is 
not  a  good  choice  for  shape  sensitivity  calculations,  due  to  the  limited  accuracy  of  the  solution. 
CTSE  requires  only  half  of  the  number  of  sample  points  required  by  CD,  thus  CTSE  only 
requires  1.5  times  more  computational  effort  than  CD.  This  coupled  with  the  fact  that  CTSE 
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doesn’t  require  changing  the  location  of  nodes  in  the  complex  plane  means  that  it  may  still  be  a 
good  choice  for  shape  sensitivity  problems. 

One  of  the  biggest  problems  with  CD  is  that  it  requires  the  user  to  change  the  location  of 
several  nodes.  Depending  on  the  method  chosen,  this  may  lead  to  elements  with  poor  aspect 
ratios,  or  complete  remeshing  of  the  domain,  or  both.  This  is  especially  true  when  the  step  size 
is  rather  large.  Since  CTSE  does  not  require  moving  the  nodes  in  the  real  plane,  remeshing  is  not 
required,  and  there  is  less  concern  over  the  aspect  ratio  of  the  elements.  Furthermore,  if  only  the 
first  order  sensitivity  is  required,  than  the  step  size  can  be  made  very  small  without  fear  of 
increasing  the  round-off  error.  This  may  be  very  useful  for  problems  in  which  external 
constraints  may  prohibit  the  domain  boundary  from  moving. 
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Tables 

Table  1.  The  Norm  of  the  Error  in  the  Stress  Solutions  for  Example  1 


Number  of  Elements 

Radial  Stress 

Tangential  Stress 

888 

5.2510e-3 

3.1852e-3 

1792 

2.6262e-3 

1.7412e-3 

17 


6184 

7.5747e-4 

5.3097e-4 

25124 

1.8608e-4 

1.3248e-4 

Table  2.  The  Computational  Time  Required  to  Solve  each  Model  for  Example  1 


Number  of  Elements 

Time  For  Solution 

888 

8.49  s 

1792 

19.48  s 

6184 

130.76  s 

25124 

2799.77  s 

Table  3.  The  Norm  of  the  Error  in  the  Sensitivity  of  the  Stress  to  the  Inner  Radius  for  Example  1 


Method 

Radial  Stress 

Tangential  Stress 

Order 

1st 

2nd 

3rd 

1st 

2nd 

3rd 

CD 

4.7265e-2 

8.1689e-2 

2.4292e-2 

8.4286e-2 

3.2213e-l 

1.4765e-2 

CTSE 

4.7268e-2 

8.1766e-2 

2.4298e-2 

8.2488e-2 

3.1587e-l 

1.4769e-2 

FD 

4.7268e-2 

8.1671e-2 

2.4295e-2 

8.7353e-2 

3.3258e-l 

1.4769e-2 

Table  4.  The  Norm  of  the  Error  in  the  First  Order  Sensitivity  of  the  Radial  Stress  as  a  Function  of  the 
Number  of  Elements  for  Example  1 _ _ _ 


Number  of  Elements 

CD 

CTSE 

FD 

Order 

1st 

2nd 

1st 

2nd 

1st 

2nd 

888 

1.2707e-l 

2.0865e-l 

1.2720e-l 

2.0856e-l 

1.2720e-l 

2.0859e-l 

1792 

8.9792e-2 

1 .53 15e-l 

8.9924e-2 

1.5214e-l 

8.9934e-2 

1.5248e-l 

6184 

4.7265e-2 

8.1689e-2 

4.7268e-2 

8.1766e-2 

4.7268e-2 

8.1671e-2 

25124 

2.2981e-2 

4.0049e-2 

2.2983e-2 

4.0088e-2 

2.2983e-2 

3.9894e-2 

Table  5.  The  Norm  of  the  Error  in  the  First  and  Second  Order  Sensitivities  of  the  Radial  Stress  as  a 
Function  of  Step  Size  for  Example  1 _ _ _ 


Step  Size 

CD 

CTSE 

FD 

Order 

1st 

2nd 

1st 

2nd 

1st 

2nd 

.0001 

4.7268e-2 

8.2959e-2 

4.7268e-2 

1.0500e-l 

4.7268e-2 

9.2942e-2 

.001 

4.7265e-2 

8.1689e-2 

4.7268e-2 

8.1766e-2 

4.7268e-2 

8.1671e-2 

.01 

4.7091e-2 

8.2920e-2 

4.7249e-2 

8.1197e-2 

4.727  le-2 

8.1748e-2 

Table  6.  The  Norm  of  the  Error  in  the  Stress  Solutions  for  Example  2 


Number  of 
Elements 

Norm  of  Error  in  Stress 
in  X 

Norm  of  Error  in 

Stress  in  Y 

Norm  of  Error  in 

Shear  Stress 

1148 

1.2615e-l 

4.3792e-2 

9.7171e-2 
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2502 

8.6330e-2 

3.2347e-2 

6.1229e-2 

8374 

4.5780e-2 

2.0857e-2 

3.9533e-2 

Table  7.  The  Computational  Time  Required  to  Solve  each  Model  for  Example  2 


Number  of  Elements 

Time  For 

1148 

10.41  s 

2502 

27.82  s 

8374 

200.00  s 

Table  8.  The  Norm  of  the  Error  in  the  First  Order  Sensitivity  of  the  Normal  Stress  in  X  as  a  Function  of 
the  Number  of  Elements  for  Example  2 _ _ _ 


Number  of  Elements 

CD 

CTSE 

FD 

Order 

1st 

2nd 

1st 

2nd 

1st 

2nd 

1148 

0.4681 

0.6378 

0.4681 

0.6379 

0.4681 

0.6378 

2502 

0.4164 

0.5673 

0.4164 

0.5675 

0.4164 

0.5674 

8374 

0.4267 

0.6572 

0.4270 

0.6581 

0.4270 

0.6577 

Table  9.  The  Norm  of  the  Error  in  the  First  and  Second  Order  Sensitivities  of  the  Radial  Stress  as  a 
Function  of  Step  Size  for  Example  2 _ _ _ 


Step  Size 

CD 

CTSE 

FD 

Order 

1st 

2nd 

1st 

2nd 

1st 

2nd 

.0001 

0.4164 

0.5679 

0.4164 

0.5726 

0.4164 

0.5688 

.001 

0.4164 

0.5673 

0.4164 

0.5675 

0.4164 

0.5674 

.01 

0.4143 

0.5609 

0.4166 

0.5707 

0.4164 

0.5674 
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Figure  1.  Convergence  of  the  Error  in  the  Radial  and  Tangential  Stress  Models  for  Example  1.  A.  The 

norm  of  the  error  in  the  radial  stress,  B.  The  norm  of  the  error  in  the  Tangential  Stress 


Radial  Stress  FEM  Solution 


Tangential  Stress  FEM  Solution 


Figure  2.  The  Numerical  Solution  of  the  Radial  and  Tangential  Stresses  for  Example  1.  A)  The  FEM  solution 
for  the  radial  stresses  B)  The  FEM  solution  for  the  tangential  stresses 
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Figure  3.  The  Error  in  the  First  Order  Sensitivity  of  the  Radial  Stress  for  Example  1.  A)  Error  in  CD,  B)  error 
in  CTSE,  C)  error  in  FD. 
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Second  Derivative  of  Radial  Stress  CD  Second  Derivative  of  Radial  Stress  CTSE 


Figure  4.  The  Error  in  the  Second  Order  Sensitivity  of  the  Radial  Stress  to  the  Inner  Radius  for  Example  1. 

A)  Error  in  CD,  B)  error  in  CTSE,  C)  error  in  FD. 
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Figure  5.  The  Finite  Element  Solution  for  the  Stresses  in  a  Disc  in  Diametrical  Compression.  A)  The  numerical 
solution  for  the  stresses  in  the  x-direction  B)  The  numerical  solution  for  the  stresses  in  the  y-direction  C)  The 
numerical  solution  for  the  shear  stress 
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Figure  6.  The  Error  in  the  Finite  Element  Method  for  Three  Different  Meshes  for  Example  2.  A)  The  error  for 
the  mesh  with  1,148  elements  and  2,357  nodes,  B)  The  error  for  the  mesh  with  2,502  elements  and  5,093  nodes,  C) 
The  error  for  the  mesh  with  8,374  elements  and  16,909  nodes. 
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1st  Derivative  of  Sigma  x  CD 
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Figure  7.  The  Error  in  the  First  Order  Sensitivity  for  Example  2.  A)  Solution  calculated  by  CD,  B)  Solution 
calculated  by  CTSE,  C)  Solution  calculated  by  FD 
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2nd  Derivative  of  Sigma  x:  CD 


2nd  Derivative  of  Sigma  x:  CTSE 
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2nd  Derivative  of  Sigma  x:  FD 
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Figure  8.  The  Error  in  the  Second  Order  Sensitivity  for  Example  2.  A)  Solution  calculated  by  CD,  B)  Solution 
calculated  by  CTSE,  C)  Solution  calculated  by  FD 
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